library(tidyverse)
library(naniar) #this is a super useful package, but slow.
library(janitor)
library(haven)
library(broom)
library(ggExtra)
library(patchwork)
library(lm.beta)
library(margins)
library(interplot)
library(modelsummary)
library(ggrepel)
library(stargazer)


#Set working directory

#setwd("c:/Documents/my/working/directory")

#skip to line 100 if you have already set up the data
#download data at https://www.europeansocialsurvey.org/data/download.html?r=8

#this takes a while (several minutes)

data_env <- read_dta("ESS8e02_2.dta/ESS8e02_2.dta") %>%
  replace_with_na_at(.vars = c("impsafe", "ipfrule", "ipbhprp", "imptrad", "trstplt", "trstplc", "trstprl", "trstlgl", "impenv", "lrscale", "ccrdprs", "trstep", "ownrdcc", "inctxff"), condition = ~.x > 11) %>%
  replace_with_na_at(.vars = c("clmchng"), condition = ~.x > 4) %>% 
  replace_with_na_at(.vars = c("wrclmch"), condition = ~.x > 5) %>% 
  mutate(auth_val = 100-(((impsafe+ipfrule+ipbhprp+ipstrgv+imptrad)-5)*4)) %>%
  mutate(auth_val=auth_val/4) %>% # Norris and Ingelhart code it as a 0-100, this just makes it a 0-25
  mutate(pop_val = (50-(trstplt+trstplc+trstprt+trstprl+trstlgl))*2) %>%
  mutate(env_imp = 6-impenv) %>%
  mutate(inctxff = 6-inctxff) %>%
  mutate(envpop = ((env_imp)/5)/((inctxff)/5)) %>%
  replace_with_na_at(.vars = c("eisced"), condition = ~.x >50) %>% #this gets rid of the category "other" in this variable, there are only 88 respondents with this value
  mutate(education = ifelse(eisced==1, "Less than Lower Secondary",
                            ifelse(eisced==2, "Lower Secondary", 
                                   ifelse(eisced==3, "Lower Tier Upper Secondary",
                                          ifelse(eisced==4, "Upper Tier Upper Secondary",
                                                 ifelse(eisced==5, "Advanced Vocational",
                                                        ifelse(eisced==6, "BA",
                                                               ifelse(eisced==7, "MA or greater", NA)))))))) %>% 
  glimpse()

# Standardize the variables of interest (not standardizing gender)

data_env <- data_env %>% 
  mutate(envpop_z = (envpop - mean(envpop, na.rm=TRUE))/sd(envpop, na.rm = TRUE)) %>% 
  mutate(auth_val_z = (auth_val-mean(auth_val, na.rm=TRUE))/sd(auth_val, na.rm=TRUE)) %>%
  mutate(lrscale_z = (lrscale - mean(lrscale, na.rm=TRUE))/sd(lrscale, na.rm=TRUE)) %>% 
  mutate(edu_z = (eisced-mean(eisced, na.rm=TRUE))/sd(eisced, na.rm=TRUE)) %>% 
  mutate(income_decile_z = (hinctnta - mean(hinctnta, na.rm= TRUE))/sd(hinctnta, na.rm=TRUE)) %>%
  mutate(rlgdgr_z = (rlgdgr-mean(rlgdgr, na.rm=TRUE))/sd(rlgdgr, na.rm=TRUE)) %>% 
  mutate(domicil_z =(domicil -mean(domicil, na.rm=TRUE))/sd(domicil, na.rm = TRUE)) %>%
  mutate(env_imp_z = (env_imp - mean(env_imp, na.rm=TRUE))/sd(env_imp, na.rm=TRUE)) %>% 
  mutate(inctxff_z = (inctxff - mean(inctxff, na.rm=TRUE))/sd(inctxff, na.rm=TRUE)) %>% 
  glimpse()

#Visually confirm that those look somewhat normalish

#hist(data_env$envpop_z) # Not particularly normal looking, but OK
#hist(data_env$auth_val_z)
#hist(data_env$lrscale_z)
#hist(data_env$edu_z)
#hist(data_env$income_decile_z)
#hist(data_env$rlgdgr_z)
#hist(data_env$domicil_z)
#hist(data_env$env_imp_z)
#hist(data_env$inctxff_z)

#create alternative dependent variable and alternative environmental collectivism variable

data_env <- data_env %>% 
  mutate(envpop_alt = (env_imp * 3)/(wrclmch+inctxff+ccrdprs)) %>%
  mutate(env_collect = wrclmch+inctxff+ccrdprs) %>% 
  glimpse()

#Visually check them

#hist(data_env$envpop_alt)
#summary(data_env$envpop_alt)
#hist(data_env$env_collect)
#summary(data_env$env_collect)

#see if and how strongly these are correlated
cor(data_env$envpop_alt, data_env$envpop, use = "complete.obs")

#it should be 0.5053127

#standardize envpop_alt and environmental collectivism
data_env <- data_env %>% 
  mutate(envpop_alt_z = (envpop_alt - mean(envpop_alt, na.rm=TRUE))/sd(envpop_alt, na.rm = TRUE)) %>%
  mutate(env_collect_z = (env_collect - mean(env_collect, na.rm=TRUE))/sd(env_collect, na.rm = TRUE)) %>% 
  glimpse()

#hist(data_env$envpop_alt_z)
#hist(data_env$env_collect_z)

saveRDS(data_env, file = "data_env_ess8.rds")
data_env <-readRDS("data_env_ess8.rds")


##create Histograms of important variables

data_env %>% 
  ggplot(aes(x=envpop))+
  geom_histogram(binwidth = .25)+
  aes(y = (..count..)/sum(..count..)) +
  scale_y_continuous()+
  theme_minimal()+
  labs(x="Environmental Populism",
       y="Proportion",
       caption = "Data From Round 8 of the ESS")

ggsave("figures/envpop_hist_in_r.png", bg="white", width = 6.5, height = 4.33)

data_env %>%
  ggplot(aes(x=(env_imp)))+
  geom_histogram(binwidth = .5)+
  aes(y = (..count..)/sum(..count..)) +
  scale_y_continuous()+
  scale_x_continuous(n.breaks = 6)+
  theme_minimal()+
  labs(x="Personal Environmental Preferences",
       y="Proportion",
       caption = "Data From Round 8 of the ESS: impenv, scale flipped")
ggsave("figures/impenv_hist.png", bg="white", width = 6.5, height = 4.33)

data_env %>% 
  ggplot(aes(x=inctxff))+
  geom_histogram(binwidth = .5)+
  aes(y = (..count..)/sum(..count..)) +
  scale_y_continuous()+
  theme_minimal()+
  labs(x="Collectivist Environmental Preferences",
       y="Proportion",
       caption = "Data From Round 8 of the ESS: inctxff, scale flipped")
ggsave("figures/inctxff.png", bg="white", width = 6.5, height = 4.33)

p1<-data_env %>%
  ggplot(aes(x=auth_val, y=env_imp))+
  geom_jitter(alpha=.5)+
  geom_smooth(method = "lm")+
  labs(x="Individual Authoritarianism",
       y="Self-Identified Environmental Importance")+
  theme_minimal()

ggMarginal(p1, type="histogram", yparams = list(binwidth = 1))

p2<-data_env %>%
  ggplot(aes(x=auth_val, y=lrscale))+
  geom_jitter(alpha=.5)+
  geom_smooth(method = "lm")+
  labs(x="Individual Authoritariansim",
       y="Political Ideology")+
  theme_minimal()

ggMarginal(p2, type="histogram", yparams = list(binwidth=1))

p3<-data_env %>%
  ggplot(aes(y=inctxff, x=lrscale))+
  geom_jitter(alpha=.5)+
  geom_smooth(method = "lm")+
  labs(y="Collectivism",
       x="Political Ideology")+
  theme_minimal()

ggMarginal(p3, type="histogram", yparams=list(binwidth=1), xparams = list(binwidth=1))



library(rstanarm)
library(tidybayes)
library(lme4)

#Run Hierarchical Models

#First the Unstandardized models

lmer_fit_unstandardized_no_inter <-lmer(envpop ~ auth_val+lrscale+eisced+gndr+domicil+hinctnta+rlgdgr+ (1|cntry/region), data=data_env, REML = FALSE)
summary(lmer_fit_unstandardized_no_inter)

lmer_fit_unstandardized <-lmer(envpop ~ auth_val*lrscale+eisced+gndr+domicil+hinctnta+rlgdgr+ (1|cntry/region), data=data_env, REML = FALSE)
summary(lmer_fit_unstandardized)

#Now the standardized models

lmer_fit_standardized_no_inter <- lmer(envpop_z ~ auth_val_z+lrscale_z+edu_z+gndr+domicil+income_decile_z+rlgdgr_z+ (1|cntry/region), data=data_env, REML = FALSE)
summary(lmer_fit_standardized_no_inter)

lmer_fit_standarized <- lmer(envpop_z ~ auth_val_z*lrscale_z+edu_z+gndr+domicil+income_decile_z+rlgdgr_z+ (1|cntry/region), data=data_env, REML = FALSE)
summary(lmer_fit_standarized)

# Standardized Models of Local Environmental Attitudes

lmer_fit_standardized_no_inter_env_imp <- lmer(env_imp_z ~ auth_val_z+lrscale_z+edu_z+gndr+domicil+income_decile_z+rlgdgr_z+ (1|cntry/region), data=data_env, REML = FALSE)
summary(lmer_fit_standardized_no_inter_env_imp)

lmer_fit_standarized_env_imp <- lmer(env_imp_z ~ auth_val_z*lrscale_z+edu_z+gndr+domicil+income_decile_z+rlgdgr_z+ (1|cntry/region), data=data_env, REML = FALSE)
summary(lmer_fit_standarized_env_imp)

# Standardized Models of Environmental Collectivism

lmer_fit_standardized_no_inter_inctxff_z <- lmer(inctxff_z ~ auth_val_z+lrscale_z+edu_z+gndr+domicil+income_decile_z+rlgdgr_z+ (1|cntry/region), data=data_env, REML = FALSE)
summary(lmer_fit_standardized_no_inter_inctxff_z)

lmer_fit_standardized_inctxff_z <- lmer(inctxff_z ~ auth_val_z*lrscale_z+edu_z+gndr+domicil+income_decile_z+rlgdgr_z+ (1|cntry/region), data=data_env, REML = FALSE)
summary(lmer_fit_standardized_inctxff_z)

#create Marginal Effects plot: Main Results

interplot(lmer_fit_standarized, var1="auth_val_z", var2 = "lrscale_z", hist =TRUE)+
  #geom_hline(yintercept = 0, linetype="dashed", col="red")+
  labs(x="Left-Right Scale (Standardized)",
       y="Marginal Effect of IA (Standardized)")+
  theme_minimal()
ggsave("figures/Standardized_Marginal_Effects.png", bg="white", width = 6.5, height = 4.33)


#create marginal effects plot: Local

interplot(lmer_fit_standarized_env_imp, var1="auth_val_z", var2 = "lrscale_z", hist =TRUE)+
  #geom_hline(yintercept = 0, linetype="dashed", col="red")+
  labs(x="Left-Right Scale (Standardized)",
       y="Marginal Effect of IA (Standardized)")+
  theme_minimal()
ggsave("figures/Standardized_Marginal_Effects_env_imp.png", bg="white", width = 6.5, height = 4.33)

#create Marginal effects plot : Collectivism


interplot(lmer_fit_standardized_inctxff_z, var1="auth_val_z", var2 = "lrscale_z", hist =TRUE)+
  #geom_hline(yintercept = 0, linetype="dashed", col="red")+
  labs(x="Left-Right Scale (Standardized)",
       y="Marginal Effect of IA (Standardized)")+
  theme_minimal()
ggsave("figures/Standardized_Marginal_Effects_inctxff_z.png", bg="white", width = 6.5, height = 4.33)


#create tables
stargazer(lmer_fit_standardized_no_inter, lmer_fit_standarized, dep.var.labels = "Environmental Populism (Standardized)",  title = "Standardized Hirerarchical Model Results", covariate.labels = c("Individual Authoritarianism", "Political Ideology", "Education", "Gender (Female)", "Rurality", "Income Decile", "Religiosity", "IA X Political Ideology"))
stargazer(lmer_fit_standardized_no_inter_env_imp, lmer_fit_standarized_env_imp, dep.var.labels = "Local Environmental Attitudes (Standardized)",  title = "Standardized Hirerarchical Model Results", covariate.labels = c("Individual Authoritarianism", "Political Ideology", "Education", "Gender (Female)", "Rurality", "Income Decile", "Religiosity", "IA X Political Ideology"))
stargazer(lmer_fit_standardized_no_inter_inctxff_z, lmer_fit_standardized_inctxff_z, dep.var.labels = "Environmental Collectivist Attitudes (Standardized)",  title = "Standardized Hirerarchical Model Results", covariate.labels = c("Individual Authoritarianism", "Political Ideology", "Education", "Gender (Female)", "Rurality", "Income Decile", "Religiosity", "IA X Political Ideology"))


#create one big table
#Note that will cause issues if your using Rtools 4.2, this website gives a great fix https://gist.github.com/alexeyknorre/b0780836f4cec04d41a863a683f91b53

stargazer(lmer_fit_standardized_no_inter, lmer_fit_standarized,lmer_fit_standardized_no_inter_env_imp, lmer_fit_standarized_env_imp,lmer_fit_standardized_no_inter_inctxff_z, lmer_fit_standardized_inctxff_z, dep.var.labels = c("Environmental Populism", "Local Environmental Attitudes", "Environmental Collectivist Attitudes"), covariate.labels = c("Individual Authoritarianism", "Political Ideology", "Education", "Gender (Female)", "Rurality", "Income Decile", "Religiosity", "IA X Political Ideology"), float.env = "sidewaystable")


# Coef Plot

modelplot(lmer_fit_standarized, coef_omit ="(Intercept)|edu_z|gndr|domicil|income_decile_z|rlgdgr_z", coef_rename = c("auth_val_z"="Individual Authoritarianism", "lrscale_z" ="Political Ideology", "auth_val_z:lrscale_z"= "IA X Political Ideology"))+
  geom_vline(xintercept = 0, linetype="dashed")+
  labs(title = "Selected Coefficent Estimates")

ggsave("modelplot.png")


#Run Stan Models (for Appendix)
# These take quite some time!

#My home computer is quadcore, so I only ran using two cores (so I could continue to use my computer for other things)
#feel free to change the cores option to anything from 1-4

#if you want to run the standardized model

#stan_fit_unstandarized <- stan_lmer(envpop ~ auth_val*lrscale+eisced+gndr+domicil+hinctnta+rlgdgr+ (1|cntry/region), data=data_env, cores=2)
#summary(stan_fit_unstandarized)



stan_fit_standardized <- stan_lmer(envpop_z ~ auth_val_z*lrscale_z+edu_z+gndr+domicil+income_decile_z+rlgdgr_z+ (1|cntry/region), data=data_env, cores=4)
summary(stan_fit_standardized)

bayesplot::color_scheme_set("brightblue")
p <- plot(stan_fit_standardized, "areas", pars = c("auth_val_z", "lrscale_z", "auth_val_z:lrscale_z", "edu_z", "gndr", "domicil", "income_decile_z", "rlgdgr_z"),
          prob = 0.5, prob_outer = 0.9)

p + ggplot2::geom_vline(xintercept = 0, linetype="dashed", col="red") +scale_y_discrete(labels = c("auth_val_z"="Individual Authoritarianism", "lrscale_z"= "Political Ideology", "auth_val_z:lrscale_z"="IA X Political Ideology", "edu_z"="Education","gndr"= "Gender (Female)", "domicil"="Rurality", "income_decile_z"="Income Decile", "rlgdgr_z"= "Religiosity"), limits = rev)


ggsave("stan_fit_standardized.png", bg="white")

#making something interesting

s <-tibble(auth_val_z=c(min(data_env$auth_val_z, na.rm=TRUE):max(data_env$auth_val_z, na.rm=TRUE)),
           lrscale_z=0,
           edu_z=0,
           gndr=0,
           domicil =3, #this is the modal category
           income_decile_z=0,
           rlgdgr_z=0,
           cntry= "DE",
           region= c("DE1"))

qi_bw <- s %>% 
  add_predicted_draws(stan_fit_standardized, ndraws = 50) %>% 
  glimpse()

s <-tibble(auth_val_z=c(min(data_env$auth_val_z, na.rm=TRUE):max(data_env$auth_val_z, na.rm=TRUE)),
           lrscale_z=0,
           edu_z=0,
           gndr=0,
           domicil =3, #this is the modal category
           income_decile_z=0,
           rlgdgr_z=0,
           cntry= "DE",
           region= c("DED"))

qi_sachsen <- s %>% 
  add_predicted_draws(stan_fit_standardized, ndraws = 50) %>% 
  glimpse()


qi <- bind_rows(qi_bw, qi_sachsen)

qi <- qi %>% 
  mutate(reg = ifelse(region=="DE1", "Baden Württemberg",
                      ifelse(region=="DED", "Sachsen", NA))) %>% 
  glimpse()


qi_sachsen %>% 
  ggplot(aes(x=auth_val_z, y=.prediction,  group = .draw))+
  geom_line()+
  theme_minimal()

qi_bw %>%
  ggplot(aes(x=auth_val_z, y=.prediction,  group = .draw))+
  geom_line()+
  theme_minimal()

qi %>%
  ggplot(aes(x=auth_val_z, y=.prediction, color = reg, group = interaction(reg, .draw)))+
  geom_line()+
  theme_minimal()

#alternate DV


stan_fit_standardized_alt <- stan_lmer(envpop_alt_z ~ auth_val_z*lrscale_z+edu_z+gndr+domicil+income_decile_z+rlgdgr_z+ (1|cntry/region), data=data_env, cores=4)
summary(stan_fit_standardized_alt)

p<- plot(stan_fit_standardized_alt, "areas", pars = c("auth_val_z", "lrscale_z", "auth_val_z:lrscale_z", "edu_z", "gndr", "domicil", "income_decile_z", "rlgdgr_z"))

p + ggplot2::geom_vline(xintercept = 0, linetype="dashed", col="red") +scale_y_discrete(labels = c("auth_val_z"="Individual Authoritarianism", "lrscale_z"= "Political Ideology", "auth_val_z:lrscale_z"="IA X Political Ideology", "edu_z"="Education","gndr"= "Gender (Female)", "domicil"="Rurality", "income_decile_z"="Income Decile", "rlgdgr_z"= "Religiosity"), limits = rev)

ggsave("stan_fit_standardized_alt.png", bg="white")



# Run hierarchical models with alternative specification of dependent variable



lmer_fit_unstandardized_no_inter_alt <-lmer(envpop_alt ~ auth_val+lrscale+eisced+gndr+domicil+hinctnta+rlgdgr+ (1|cntry/region), data=data_env, REML = FALSE)
summary(lmer_fit_unstandardized_no_inter_alt)

lmer_fit_unstandardized_alt <-lmer(envpop_alt ~ auth_val*lrscale+eisced+gndr+domicil+hinctnta+rlgdgr+ (1|cntry/region), data=data_env, REML = FALSE)
summary(lmer_fit_unstandardized_alt)

interplot(lmer_fit_unstandardized_alt, var1="auth_val", var2 = "lrscale", hist =TRUE)+
  #geom_hline(yintercept = 0, linetype="dashed", col="red")+
  labs(x="Left-Right Scale",
       y="Marginal Effect of IA")+
  theme_minimal()
ggsave("figures/Unstandardized_Marginal_Effects_alt.png", bg="white", width = 6.5, height = 4.33)



#Now the standardized models of Environmental Populism (alt)

lmer_fit_standardized_no_inter_alt <- lmer(envpop_alt_z ~ auth_val_z+lrscale_z+edu_z+gndr+domicil+income_decile_z+rlgdgr_z+ (1|cntry/region), data=data_env, REML = FALSE)
summary(lmer_fit_standardized_no_inter_alt)

lmer_fit_standarized_alt <- lmer(envpop_alt_z ~ auth_val_z*lrscale_z+edu_z+gndr+domicil+income_decile_z+rlgdgr_z+ (1|cntry/region), data=data_env, REML = FALSE)
summary(lmer_fit_standarized_alt)

interplot(lmer_fit_standarized_alt, var1="auth_val_z", var2 = "lrscale_z", hist =TRUE)+
  #geom_hline(yintercept = 0, linetype="dashed", col="red")+
  labs(x="Left-Right Scale (Standardized)",
       y="Marginal Effect of IA (Standardized)")+
  theme_minimal()
ggsave("figures/Standardized_Marginal_Effects_alt.png", bg="white", width = 6.5, height = 4.33)

#Standardized Models of Environmental Collectivism (alt)

lmer_fit_standardized_no_inter_env_collect <- lmer(env_collect_z ~ auth_val_z+lrscale_z+edu_z+gndr+domicil+income_decile_z+rlgdgr_z+ (1|cntry/region), data=data_env, REML = FALSE)
summary(lmer_fit_standardized_no_inter_env_collect)

lmer_fit_standardized_env_collect <- lmer(env_collect_z ~ auth_val_z*lrscale_z+edu_z+gndr+domicil+income_decile_z+rlgdgr_z+ (1|cntry/region), data=data_env, REML = FALSE)
summary(lmer_fit_standardized_env_collect)

# Coef Plot

modelplot(lmer_fit_standarized_alt, coef_omit ="(Intercept)", coef_rename = cc("auth_val_z"="Individual Authoritarianism", "lrscale_z"= "Political Ideology", "auth_val_z:lrscale_z"="IA X Political Ideology", "edu_z"="Education","gndr"= "Gender (Female)", "domicil"="Rurality", "income_decile_z"="Income Decile", "rlgdgr_z"= "Religiosity"))+
  geom_vline(xintercept = 0, linetype="dashed")+
  labs(title = "Selected Coefficent Estimates")

ggsave("modelplot_alt.png")

#create table

stargazer(lmer_fit_standardized_no_inter_alt, lmer_fit_standarized_alt,lmer_fit_standardized_no_inter_env_imp, lmer_fit_standarized_env_imp,lmer_fit_standardized_no_inter_env_collect, lmer_fit_standardized_env_collect, dep.var.labels = c("Environmental Populism", "Local Environmental Attitudes", "Environmental Collectivist Attitudes"), covariate.labels = c("Individual Authoritarianism", "Political Ideology", "Education", "Gender (Female)", "Rurality", "Income Decile", "Religiosity", "IA X Political Ideology"), float.env = "sidewaystable")